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ABSTRACT 

The objective of the present work was the structural optimization of thin shell structures 
that are subjected to stress and displacement constraints. In order to accomplish this, the 
structural optimization computer program DESAPl was modified and improved. In the static 
analysis part of the DESAPl computer program the torsional spring elements, which are used to 
analyze thin, shallow shell structures, were eliminated by modifying the membrane stiffness 
matrix of the triangular elements in the local coordinate system and adding a fictitious 
rotational stiffness matrix. This simplified the DESAPl program input, improved the accuracy 
of the analysis, and saved computation time. In the optimization part of the DESAPl program 
the stress ratio formula, which redesigns the thickness of each finite element of the structure, 
was solved by an analytical method. This scheme replaced the iterative solution that was 
previously used in the DESAPl program, thus increasing the accuracy and speed of the design. 

The modified program was used to design a thin, cylindrical shell structure with optimum 
weight, and the results are reported in this paper. 



1. INTRODUCTION 


1.1 General 

In the past, discrete optimality criteria have been derived for a number of design conditions 
including strength, static stiffness, dynamic stiffness, static stability, and aeroelastic constraints. 
The computer programs DESAP [1], DESAP2 [2], and FASTOP [3], as well as OPSTAT, 
OPTCOMP, OPTIM, ASOP [4], and others, use the discrete optimality criteria approaches as a 
basis for structural optimization. These programs are based on the finite element method of 
analysis and can optimize isotropic, anisotropic, and layered composite structures. The codes 
can handle over 1000 design variables and a comparable number of analysis variables, which are 

the degrees of freedom. 

This paper examines the structural optimization of thin, shallow shells by using the 
DESAP 1 structural optimization program. DESAP1 was developed to design structures with 
linear elastic material behavior. The total weight of the structure is minimized by computing 
the element sizes, that is, the cross-sectional areas for trusses and beams and the thicknesses for 
plates and shells, under certain constraints. The static analysis of the design is computed by 
using the finite element method. For purpose of analysis the SAPIV finite element program [5] 
was modified for use in the DESAP1 computer program. The synthesis algorithm of DESAP1 
consists of iterative processes. Each iterative process comprises three steps: (1) static analysis 
of the current design, (2) comparison and evaluation of the results of the static analysis, and 
(3) redesign of the structure by using information from the previous two steps. 

Two kinds of constraints are used in DESAP1: primary constraints and secondary 
constraints. The primary constraints are displacements and stresses with upper limits. The 
stress constraints are failure criteria, or local instability criteria, or both. If stress constraints 
are used in the redesign of the structure, the stress ratio method [1, 6, 7] is applied to drive the 
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final design to the fully stressed design. In the fully stressed design each element is assumed to 
reach the allowable stress under at least one load condition. The fully stressed design always 
coincides with the minimum-weight design for statically determined structures, but the design 
need not be fully stressed at optimum for statically indeterminate structures, in general. For 
displacement constraints the redesign procedure [6] is based on an optimality criteria method, 
and the element sizes (thicknesses) are obtained for the optimum weight design. 

The secondary constraints consist of the minimum allowable element sizes and the size 
proportional constraints, whereas the element sizes have a prescribed ratio. In the DESAP 1 
program it is assumed that the static loading is independent of the element sizes. Also, during 
the design procedure the layout of the structure is not changed. 

This study involved the following steps: 

(1) In the static analysis part of the DESAPl program the torsional spring elements, which 
are used to analyze thin, shallow shell structures, were eliminated by modifying the membrane 
stiffness matrix of the triangular elements in the local coordinate system and adding a fictitious 
rotational stiffness matrix as suggested by Zienkiewicz [8, 9]. 

(2) In the optimization part of the DESAPl computer program the iterative solution of a 
fourth-order equation, which redesigns the thickness of each finite element of the structure, was 
replaced by an analytical solution. 

(3) A thin shell structure was designed to minimize its weight and was subjected to stress 
and displacement constraints. 

1.2 Element properties 

The weight Wj of an element i can be written as 

w i = Pi tj i = 1,..,I 
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where pj is the unit weight and tj is the size (thickness) of the element i. Because it is not 
desirable to have each tj as an independent variable, the form 

H = n i D m 

is introduced for i = and m = < I, where D m are independent variables and n ; is 

the design variable fraction of the element. For each element, nj and m must be defined [l]. 
The stiffness matrix [Kj] of element i can be written as 

w = [KfVi + (Kfvr 

where [k[^] is the unit stiffness matrix due to the action of the direct stresses, [K: J ] is the unit 
stiffness matrix due to the bending or torsion, and nj is the inertia exponent, which is 
determined by the relationship between the size and moment of inertia of the element as follows: 

Ii = J,tf 

where Jj is the unit moment of inertia. 

The vector of the internal forces {Nj} of node i is computed from the previously computed 
element nodal displacements {uj} as follows: 

(N,) = [S,]{Uj} + {T,} 

where [S-] is the recovery matrix and {Tj} is the force vector [1]. 


1.3 Stress-ratio formula 

If stress constraints are imposed in the DESAP 1 program, the stress ratio method is used. 

In order to obtain the redesign stress ratio formula, the Von Mises yield criterion is applied at 
each element of the plate or shell structure. 

First, the assumption is made that a single load condition is imposed in the structure. The 
yield criterion f for an element i has the general form: 

f({Ni>, {N*}, t) = 1 *' 
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where {NJ is the vector of the internal forces on the nodes of an element i due to a single load 

condition; {NJ is the vector of the allowable forces; and t is the thickness of an element i of 
the plate or shell structure. 

If the improved design is to be fully stressed and {Nf} and t' are the improved values, 
the following form must be satisfied: 


f ({N/}, {NJ, t') = 1 (2) 

If {Nj} is known, the improved value of the thickness t' is calculated from Eq. (2). The 

calculation of {N/} is obtained by inverting the structural stiffness matrix. In the DESAPl 

program this inversion is impractical for large structures because the banded form of the 

structural stiffness matrix is lost after its inversion. For this reason it is assumed that the nodal 

forces {N/} = {NJ do not change. Substituting the values of {N/} into Eq. (2) gives the 
following equation: 


f«N,>, {N*}, f) = i 

The solution procedure consists of solving Eq. (3) for f load conditions at a time and 


( 3 ) 


then choosing the largest value for the improved thickness. This procedure is the stress ratio 
method of redesign as it is used in the DESAPl program. Equation (3) is an approximation and 
must be applied iteratively, updating {N,} each time, before a full, stressed design is obtained. 

In order to check the design process, the following two controls were established: 

(1) The design is critical if 


( l ~ d ) < Rmax < (1 + d) 

where d is a small parameter described by the user of the program that expresses the desired 
width of the critical design band; R m „ is equal to max (E&/D„) , where D m and D^, are 
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,h« previous end the new values, respectively, and % is the largest value of the following two 
equations: 

Di = t'/ni 

where nj is a parameter described by the user for an element i, or 

££ 

where is the prescribed lower thickness of the secondary constraints. When equal or 
proportional sire constraints are used, is defined by choosing the largest value of the 
following two equations: 

D' = max (t'/nO 


or 


%= D m 

(2) The design is fully stressed if it satisfies the foUowing two conditions simultaneously: 
(1 - d) < R max < (1 + d) (critical design) 


and 


(1 - b) < Rnan < <> + d > 
where R m ;„ is equal to min (D m D m ). 


1.4 Hencky-Von Mises failure criterion 

In this study the plate or shell structure was assumed to consist of material that is isotropic 
and homogeneous. The Hencky-Von Mises failure criterion was used to obtain the redesign 

stress ratio formula. The failure function f of Eq. (1) takes the following form: 

* * * * , * ( 4 ) 

f = (S x /S x ) 2 + (Sy/Sy) 2 - S X S y /S x S y + (S X y/S X y) - 
where S x and S y are the normal stresses at a point in the element, S xy is the shear stress 
point in the element, and S* x , S y , and S* xy are the allowable normal and shear stresses. 
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The basic assumption was that the resulting internal forces remain unchanged during the 


redesign: 

N' = N X 
n; = N y 
Nxy = N xy 

M; = My 

Ky = Mxy 

where N' and M' are the values of the improved design and N and M are the values of 
the previous design. 

The internal stresses that are developed in the shell under the external static load consist of 


the membrane components and the bending components. The membrane stresses 
S* S* and S * at a point are given by 


s; = vt' 

(5a) 

s; = N y /t* 

(5b) 

s; y = Nvt' 

(5c) 

where S' is the value of the improved design and t' is the value of the 
an element of the shell structure. The allowable membrane stresses are 

s; = N> 

S‘ y = N*/t 
S'x, = N* xy /f 

The bending stresses S x , Sy, and S xy at a point are given by [10] 

improved thickness of 

s; = ieiv^/t' 2 

(6a) 

s; = ±6My/t' 2 

(6b) 
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The allowable bending stresses are 


K, = A ' 2 


(6c) 


S* x = 6M*/t J 
S* y = 6M;/t 2 


S* = 6M* /t 2 

xy xy' 


Therefore, 

in a shell structure the stresses at a point are given by the following forms: 



” (®x/®x)membrane (®x/®x)bending 

(7a) 


(®y/®y) (®y/®y)membrane (®y/®y) bending 

(7b) 


/ * / * / * 
(®xy/®xy) (®xy/®xy) membrane (®xy/®xy)bending 

(7c) 

From Eq. (7a) 


(s;/s* ) = (N,/N j (t/f) ± (M,/M>/f) 2 

(8a) 

Similarly, 


(Sy/Sj ) = (Ny/Ny) (t/f) i (M,/M y ) ( t/ 1 ' ) ^ 

(8b) 

and 


(S^r/Sjy) = (N X y/N^) (t/f) ± (Myy/M^Kt/f) 5 

(8c) 


Multiplying Eqs. (8a) and (8b) gives 

(s;/s;)/(s>;> = KVN*Ht/f) ± (»^/iv<) (t/t-) 2 ] (8 

x |(N y /N*)(t/t') ± (My/Mj)(t/l') 2 ] 

Substituting the values S x /S*, Sy/S*, and S^ y /S xy from Eqs. (8) into Eq. (4) yields the 
following equation: 
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l(N x /N>/t') ± (M,/J<)(t/t') 2 ] 2 + |(N y /N;)(t/t') ± (M y /M y , )(t/l') 2 ] J 
+ ICV'OW) ± (K,/K> (t/t') 2 ] 1 - [(N X /N>A') ± (M^Mt/l') 2 ] 

x KVNjHt/t') ± (Mj/My'Ht/f) 2 ] = 1 

Simplifying this equation gives 

(t'/t) 4 - C„(t'/t) 2 =F C x (t'/t) -C = 0 (9) 

where the minus sign preceding C x is applicable to the upper surface of the plate or shell, and 
the plus sign, to the lower surface. In Eq. (9) the following notation was used: 

C» = (N x /Nx‘) 2 + + (N.X/ - N x N y /N*N y (10a) 

C x = 2 [(N x /N‘)(M x /M;) + (N y /N;)(M,,/M;) + (N*/N*) 

- [(Ny/N.lHMy/l^) + (N y /N y ) (Mxy/M^y)] 

and 

C = (My/I^) 2 + (My/I^) 2 + (Myy/M^y) 2 ~ 1^/1^) ( ' 0C) 


Further investigation into Eqs. (10a) and (10c) (see Appendix) gives >. 0 and C 0. 

These inequalities were used to find the positive and real solution of the quartic Eq. (9). 

In order to solve the quartic equation 

X 4 - C M X 2 =F C x X - C = 0 ( 11 ) 

where X = t'/t, it was reduced to the resolvent cubic equation [11]. The resolvent cubic 
equation of the quartic is 

y 3 + Ay 2 + By + D = 0 ( 12 ) 
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where A = -C^, B = 4C, and D = -C* - 40^0. The analytical solution of the roots of the 
cubic Eq. (12) according to Cardan’s method is described by Borofsky [12]. By the analytical 
solution of the quartic equation, a much greater accuracy and computational efficiency was 
achieved than when an iteration method was used to obtain the solution. 

2. SHELL STRUCTURES 

The original DESAP 1 computer program uses torsional spring elements, which are springs 
that are defined as normal to the flat shell surface, to eliminate the singularity of the total 
stiffness matrix in the global coordinate system. Because it is time consuming to define the 
normal direction at each point on the shell structure surface, the DESAP 1 subroutines were 
modified so that thin shells could be designed without using the torsional spring elements. In 
order to explain the procedure that was followed, shell structures and finite element theory are 
reviewed briefly. 

The classical theory of shell structures is discussed by Flugge [13]. When applying the finite 
element method to shell problems, it is assumed that the behavior of a continuously curved 
surface can be represented by the behavior of a surface that is built up with small flat elements. 
In any shell structure the elements generally will be subjected to both bending and in-plane 
forces [14]. For example, consider typical triangular flat elements that are subjected 
simultaneously to in-plane and bending actions (Figs. 1 and 2). In Fig. 1, u n and v n (n = i, j, k) 
are the nodal displacements, and U n and V n are the corresponding nodal forces. In Fig. 2, 
W n , ^xn> ^yn ( n = L j >k) are the nodal displacements, and W n , M^, and My n are the 
corresponding nodal forces. Taking first the in-plane (plane stress) action, the state of the strain 
is uniquely described in terms of the Uj and Vj displacement of each typical node i in the x' 
and y' local coordinate directions, respectively. The minimization of the total potential 
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energy led to the stiffness matrices described by Zienkiewicz [8], which relate nodal forces F p to 
displacement parameters 

F p = K p d p 

For a triangular element K p is a 6-by-6 matrix [7]. 

Similarly, when bending is considered, the state of strain is given uniquely by the nodal 
displacement w* and the two rotations 6 ^ and at each node i (Fig. 3). The 
corresponding nodal forces are given by 

F b = K b d b 

For a triangular element K b is a 9-by-9 matrix [7]. 

Before combining the in-plane and bending stiffnesses note (1) that the displacements 
prescribed for in-plane forces do not affect the bending deformations and vice versa and (2) that 
rotation 0 zi about the z' axis is not involved in the deformations in either node. Combining 
the membrane and bending actions and introducing 0 zi and its associate couple M zi gives for 
an element node i the following nodal displacements dj and nodal forces Fj: 

d* = lu- v- w* 0 • 0 * $ 
ph i i xi v yi v zi| 

and 

F t =ju i V i W i M xl M ri M,J T 

For an element 

F = kd ( 13 ) 

where for a triangular element containing nodes i, j, and k the element force and displacement 
vectors are 

F = h Fj F k | T 
d = |d 1 d,d k l 
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and the element stiffness matrix in the local coordinates is 

K = K p + K b ( 14 ) 

The total stiffness matrix K is an 18-by-18 matrix and is illustrated in Fig. 3. 

2.1 Replacement of torsional spring elements 

In a shallow shell structure, if all the elements meeting at a node i are coplanar, numerical 
problems due to the singularity of the stiffness matrix arise in the DESAP 1 program. Because it 
was assumed that the moment and the stiffness are zero in the 0 2 - direction in the local 
system (Fig. 1), the sixth row and column of each submatrix of the total stiffness matrix 
contain zeroes [7]. If the set of all equilibrium equations is considered at the point i, in the local 
system, six equations result, of which the sixth equation is 

M zi = Oui + Ovi + Owi + 00^ + O0 yi + O0 zi (15) 

or 0 — 0. The difficulty persists when the six equilibrium equations at the point i are 
transformed to global coordinates. Because these equations will still have a singular matrix, a 
solution cannot be obtained [8, 9, 15]. The following three techniques for eliminating the 
singularity of the stiffness matrix have been used in existing literature: 

(1) References 8,15, and 16 eliminate the sixth row and column from the stiffness matrix to 
obtain a nonsingular 5-by-5 stiffness matrix. However, it is impractical to apply this technique 
to DESAP 1 because programming difficulties are encountered that would require an extensive 
revision of the DESAP 1 source program. 

(2) The DESAP 1 and SAPIV computer programs use the torsional spring elements to 
eliminate the singularity of the stiffness matrix. The spring element has a torsional stiffness 
that is normal to the shell surface (Fig. 4). Two methods are used to specify the direction of the 
spring elements in DESAPl (Figs. 5 and 6). In the first method the direction is determined by 
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the structural node N and a second node I (Fig. 5). Node I may be a structural node or a special 
node. In the latter case the degrees of freedom of node I should be suppressed on the node. In 
the second method the direction of the element is taken as perpendicular to the lines IJ and KL 
(Fig. 6). The points I, J, K, and L may be structural nodes or special points (with suppressed 
degrees of freedom). The torsional spring elements give rise to additional terms along the 
diagonal of the element stiffness matrix K [7], These new terms have values that are equal to 
the rotational stiffness of the spring. Transforming the matrix K of each element to the global 
system and adding all the matrices give a nonsingular stiffness matrix of the structure. 

Although the use of torsional spring elements in eliminating the singularity is feasible, it is 
both difficult and time consuming to find the directions that are normal to the shell structure at 
the nodal points in order to define the new nodes and the torsional stiffnesses associated with 
them. Therefore, discovering a simpler way to overcome the use of boundary elements would 
save time and make the program easier to use. 

(3) References 8, 9, and 17 eliminate the singularity of the stiffness matrix by adding a 
fictitious rotational stiffness matrix in the membrane stiffness matrix K p . Here it is assumed 
that a nodal rotation 6 zl about the % f local axis of any one triangular element node is 
responsible only for the development of resisting couples M 2 j, M 2 j, and at the three 
element nodes and does not produce any other reactions. In order to ensure static equilibrium, 
the sum of the couples M zi , M z j, and M zk is always taken to be zero. Zienkiewicz et al. [8, 9] 
determined that satisfactory results are obtained if these couples are proportional to the modulus 
of elasticity E and the volume At of the triangular element, where A is the area and t is 
the thickness of the element. They suggested the following moment-rotation relationship: 
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1 -0.5 -0.5 



Hj 

= aEAt 

-0.5 1 -0.5 

• 




-0.5 -0.5 1 


•a 


( 16 ) 


where a is an undetermined coefficient. A realistic value of a was estimated to be 
approximately 0.03 [9, 17]. A displacement error of about 10 percent is caused by introducing a 
value as large as 1.0. The displacements for very small values of a are nearly exact. However, 
for practical purposes, extremely small values of a are possible only when a large 
computational precision is available. In the present research a value of a = 0.03 was used so 
that Eq. (16) can be written as 


Hi 


B -B/2 -B/2 


Ki 


= 

-B/2 B -B/2 

• 


Mzk 


-B/2 -B/2 B 


0*k 


(17) 


where B = 0.03EAt. Adding the stiffness coefficients B and -B/2 of Eq. (17) in the 
appropriate positions in the membrane stiffness matrix K p in the local system yields a new 
membrane stiffness matrix of the element [7]. The combination of the updated membrane 
stiffness matrix find the bending stiffness matrix gives the total stiffness matrix in the local 
system for a triangular element i, j, k (Fig. 3). (In Fig. 3, C = B/2.) Transforming the total 
stiffness matrix from the local to the global system of each element and adding them gives a 
nonsingular structural stiffness matrix. This procedure eliminates the singularity of the 
structural stiffness matrix in the global system for shallow shell structures in the DESAP 1 
computer program. 
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3. EXAMPLE PROBLEM: THIN SHELL STRUCTURE 


Consider a thin, cylindrical shell structure with an angle of 10° (Figs. 7 and 8) that is 
subjected to a concentrated load of 128 kips at the center. The shell, 36 ft long and 30 ft wide, 
is simply supported at the two opposite edges and is free at the other two edges. The material 
behavior of the structure is linear elastic, isotropic, and homogeneous. The design of the shell 
structure was obtained by the modified DESAP 1 program. Because of the symmetry of the shell 
structure, only one-fourth of the shell (Fig. 9) is required to employ the finite element method. 

Young’s modulus of elasticity of the material is 4.38xl0 8 lb/ft 2 , the specific weight is 
360 lb/ft 3 , and the Poisson’s ratio is 0.3. The concentrated load on the quarter shell structure is 
32 kips. The design commenced with a uniform thickness t for all the elements equal to 
1.14 ft. The minimum size constraint for all the element thicknesses was 0.1 ft. All the plate 
elements were sized independently. 

Both displacement and stress constraints were used. Only the displacement constraints 
were used first to allow for the thickness of the element to converge faster; they were applied at 
the center of the shell structure. The magnitude of the displacement constraint is 0.055 ft in the 
z direction. The stress constraints are added later, and the magnitude of the allowable stress in 
tension is 25 000 lb/ft 2 . The allowable stress in compression is also 25 000 lb/ft 2 . 

The design is acceptable if it meets two criteria [1, 6]. The first criterion is satisfied (1) if 
the design is fully stressed, as discussed earlier, and (2) if the displacement constraints are not 
violated; that is, Q max <_ (1 + d), where Q max is the displacement ratio, which is defined as the 
displacement of a node over the maximum allowable displacement at the node, and d is a small 
parameter described previously. The second criterion is satisfied (1) if the design is displacement 
critical (i.e., (l — d) < Qmax < (1 -f d), (2) if the stress constraints are not violated 
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(i.e., R max < (1 + d), where R max is defined as before), and (3) if the optimality criterion is 
satisfied for the displacement constraints. Figures 10 and 11 show the results of the computed 
output. 

Figure 10 shows the weight of the structure versus the number of critical designs. Note that 
the weight decreases under the displacement constraints. When after 42 cycles the stress 
constraints are added, the structural weight increases suddenly and then begins to decrease again 
until, at 44 cycles, the design becomes optimal. Figure 11 shows the thickness of selected 
elements versus the number of critical designs. In the final stages, when the design approaches 
the acceptable criteria, the thickness of the elements decreases, whereas at the intermediate 
stages it may increase or decrease. Figure 11 shows that, as expected, the elements closest to 
the central load elements have a greater thickness than the more distant elements that have 
achieved the minimum allowable thickness. The larger stresses are developed close to the central 
load elements. 


4. CONCLUSIONS 

The example problem demonstrates that optimal design of thin shells of arbitrary shape can 
be accomplished by a general-purpose synthesis program such as DESAP1. The major deficiency 
of the present algorithm is that it requires a large number of design cycles to reach convergence 
to the final design, which is typical of other structural optimization algorithms. The economy of 
computation could be improved considerably by using a uniform scaling operation. 

Because the displacement-constraint design is based on the exact optimal criterion and thus 
has better convergence characteristics, it is important in running the program to use the 
displacement constraints first and then impose the stress constraints. Once the design has 
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converged under the displacement constraints, adding the stress constraints will result in only a 
few more design cycles. 

If the design has stress constraints only, it is still worthwhile to impose contrived 
displacement constraints on the initial stages of the design procedure. These constraints should 
be replaced by the stress constraints after the initial design has converged. 
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APPENDIX — INVESTIGATION OF COEFFICIENTS C H and C 
OF REDESIGN STRESS RATIO FORMULA 
The stress ratio formula is a fourth-order equation in the variable X = t'/t, where t is 
the current thickness at the element and t' represents the improved thickness 

X 4 - C^X 2 tC x X-C = 0 (11) 

where the coefficient has the form 

C « = ( N x/N *) 2 + (Ny/Ny*) 2 + (N^/N^y) 2 - N x N y /N*N* (1 ° a) 

where N x , N y , and N xy are the normal and shear forces and N x , N*, and N xy are the 
allowable normal and shear forces. Equation (10a) can be written as 

C xx = a 2 + b 2 + c 2 - ab 

where a = N x /N*, b = N y /N*, and c = N xy /N*y. If ab is a nonpositive number (i.e., 
a b <_0), then is always a positive number. If ab is a nonnegative number (i.e., ab > 0), 

then 


-ab > -2ab 


Adding a 2 + b 2 + c 2 to both sides of this inequality gives 

a 2 + b 2 + c c - ab > a 2 + b 2 + c 2 -2ab 


or 


a 2 + b 2 + c 2 - ab > (a - b)* + c : 


2 , _2 


Therefore, 


a 2 -(- b 2 + c 2 - ab > 0 


Consequently, C n is always a nonnegative real number if ab > 0. 
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The coefficient C of Eq. (11) has the form 


C = + (My/M,') 2 + (Mjiy/Mfy) 2 - (My/M^HMy/M;) 

where M x , My, and M xy are the bending and twisting moments and M^, and M^ y are the 
allowable bending and twisting moments as described previously. Again it can be shown that 
the coefficient C is always a nonnegative real number. 

This reasoning was used for the analytical solution of Eq. (11). 


19 



REFERENCES 


1. J. Kiusalaas, and G.B. Reddy, DESAPl— A Structural Design Program with Stress and 

Displacement Constraints. Vols. 1-3. NASA CR-2794 to 2796 (1977). 

2. J. Kiusalaas and G.B. Reddy, DESAP2— A Structural Design Program with Stress and 

Bending Constraints. Vols. 1-3. NASA CR-2797 to 2799 (1977). 

3. K. Wilkinson, et al., FASTOP: A flutter and strength optimization program for lifting- 

surface structures. J. Aircraft , 14, 581-587 (1977). 

4. V.P. Venkayya, Structural Optimization: A review and some recommendations. Int. J. 

Numer. Methods Eng., 13, 203-228 (1972). 

5. K.J. Bathe, E.L. Wilson, and F.E. Peterson, SAPIV: A Structrual Analysis Program for 

Static and Dynamic Response of Linear Systems. University of California, Berkeley, CA 
(1974). 

6. J. Kiusalaas, Minimum Weight Design of Structures via Optimality Criteria. NASA TN 

D-7115, 31-41 (1972). 

7. P.K. Gotsis, Optimization of Thin Shell Structures by the Finite Element Method. M.S. 

Thesis. Engineering Science and Mechanics Department, Pennsylvania State University, 
University Park, PA (1982). 

8. O.C. Zienkiewicz, The Finite Element Method . Third ed. McGraw-Hill Book Company, 

New York, 335-337 (1977). 

9. O.C. Zienkiewicz, C. Parekn, and I.P. King, Arch Dams Analysed by a Linear Finite Element 

Shell Solution. Proc. Symp. Arch Dams, Inst. Civil Eng., London, UK, 19-36 (1968). 

10. S. Timosenko and S. Woinowsky-Krieger, Theory of Plates and Shells . McGraw-Hill Book 

Company, New York, 37-52 (1959). 


20 


11- L.W. Griffiths, Introduction to The Theory of Equations . John Wiley & Sons, Inc., 

New York, 33-42 (1947). 

12. S. Borofsky, Elementary Theory of Equations. The Macmillan Company, New York, 

115-130 (1961). 

13. W. Flugge, Stresses in Shells. Springer-Verlag, New York (1973). 

14. I. Holland, Fundamentals of the finite element method. Comput. and Struct., 4, 3-15 (1974). 

15. R.W. Clough and E.L. Wilson, Dynamic finite element analysis of arbitrary thin shells. 

Int. J . Solids Struct ., 7, 33-56 (1971). 

16. R.W. Clough and C.P. Johnson, A finite element approximation for the analysis of thin 

shells. Int . J. Solids Struct ., 4, 43-60 (1968). 

17. D. G. Aswell and R. H. Gallagher, Finite Elements for Thin Shells and Curved Members. 

John Wiley & Sons, New York, 245-263 (1974). 




Figure 1. — ln-plane forces and deformations in local 
system X' Y* Z'. 


Figure 2. — Bending forces and deformations in local 
system X 1 Y' Z*. 
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Figure 3. — Stiffness matrix K and fictitious stiffness coefficients. 



Figure 4. — Direction of boundary elements in typical 
shell structure. 



Figure 5. — First method of determining 
boundary elements. 
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Figure 10. — Optimum design of thin, fiat shell 
structure. 
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Figure 1 1 .—Thickness of element versus number of 
designs. 
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